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ABSTRACT  (UNCLASSIFIED) 


The  tracking  problem  of  unknown  marine  platforms  using  passive  sonar 
measurements  is  generally  referred  to  as  Target  Motion  Analysis  (TMA) . 
This  report  describes  the  application  of  a  Maximum  Likelihood  Estimation 
(MLE)  method  to  obtain  position  and  velocity  estimates  of  a  platform 
when  bearing  angle  and  frequency  measurements  are  available  from  a 
passive  sonar  system.  The  frequency  measurements  are  related  to  one  or 
more  cardinal  frequency  peaks  from  the  radiated  frequency  spectrum  by 
the  platform. 

The  MLE  method  offers  the  opportunity  to  use  a  multi-leg  model,  i.e.  the 
platform  to  be  localized  is  assumed  to  move  according  to  a  piecewise 
linear  track,  where  each  part  is  referred  to  as  a  leg.  On  each  leg 
constant  course  and  speed  is  assumed.  Additionally,  bearing  and 
frequency  measurements  related  to  bottom  reflections  of  the  acoustical 
signals  can  be  used. 

In  ref.  [Gmellg  Meyling,  1989-1)  a  Newton- type  optimization  method  using 


first  and  second  order  derivative  information  of  the  residual  functions  □ 
is  proposed.  By  using  the  analytic  expressions  of  the  derivatives  as 

decrlbed  in  this  report  a  major  reduction  of  computation  time  is  - 

accomplished. 
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SAMENVATT1NG  (ONGERUBRICEERD) 

In  het  algemeen  wordt  het  volgprobleem  m.b.t.  onbekende  platforms  op  zee 
waarbij  passieve  sonar  metingen  worden  gebruikt  Target  Motion  Analysis 
(TMA)  genoemd.  Dit  rapport  beschrijft  de  toepassing  van  een  Maximum 
Likelihood  schattlngsmethode  (MLE)  om  de  bewegingsparameters  van  een 
platform  te  schatten  uit  peilings-  en  frequentiemetingen  die  m.b.v.  een 
passief  sonarsysteem  worden  bepaald.  De  frequentiemetingen  zijn 
gerelateerd  aan  een  of  meer  pieken  uit  het  door  het  platform  uitgezonden 
frequentiespectrum. 

De  MLE  methode  geeft  de  mogelijkheid  een  multi-leg  model  te  gebruiken, 
d.w.z.  het  platform  wordt  verondersteld  te  varen  volgens  een  stuksgewijs 
lineaire  baan,  waarbij  elk  stuk  een  poot  genoemd  wordt.  Op  elke  poot 
wordt  een  constante  koers  en  vaart  verondersteld.  Tevens  is  het  mogelijk 
metingen  te  gebruiken  die  afkomstig  zijn  van  bodemgereflecteerde 
akoestische  signalen. 

In  ref.  [Gmelig  Meyling,  1989-1]  is  een  Newton-type 

optimalizatiemethode  voorgesteld  die  gebaseerd  ij  op  de  eerste  en  tweede 
afgelelden  van  de  resldufuncties .  Een  belangrijke  reduktie  in  de 
benodigde  rekentijd  wordt  bereikt  door  gebruik  te  maken  van  de  in  dit 
rapport  beachreven  analytische  uitdrukkingen  voor  de  afgeleiden. 
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1  INTRODUCTION 

The  Cracking  problem  of  unknown  marine  platforms  using  sonar 
measurements  Is  generally  referred  to  as  Target  Motion  Analysis  (TMA) . 
This  report  describes  the  application  of  a  Maximum  Likelihood  Estimation 
(KLE)  method  to  obtain  position  and  velocity  estimates  of  a  platform 
when  bearing  angle  and  frequency  measurements  are  available  from  a 
passive  sonar  system,  such  as  a  towed  array.  The  frequency  measurements 
are  related  to  one  or  more  cardinal  frequency  peaks  from  the  radiated 
frequency  spectrum  by  the  platform. 

The  estimation  problem  Is  formulated  as  a  Maximum  Likelihood  Estimation 
problem  (MLE)  which  is  both  non-linear  and  poorly  conditioned.  In 
literature  [Lindgren  1978,  Aidala  1979,  1982,  1983]  much  attention  is 
paid  to  the  application  of  Kalman  Filters  (KF)  to  solve  bearings-only 
TMA  problems.  In  [Nardone,  1984]  the  fundamental  properties  if  the 
bearins-only  TMA  problem  are  discussed.  An  extended  Kalman  filter  for 
the  bearing-  and  frequency  measurement  TMA  problem  is  proposed  in 
[Ockeloen  and  Hlllemsen,  1982].  However,  in  many  situations  Kalman 
filters  suffer  from  unacceptable  bias  and  slow  convergence  caused  by 
inappropriate  linearization  of  the  non-linear  TMA  equations. 

By  using  a  proper  numeric  optimization  method  to  obtain  an  ML  estimate 
the  disadvantages  of  Kalman  filters  have  been  overcome  at  the  cost  of 
more  computational  effort. 

The  optimization  method  in  [Gmelig  Meyling,  1989]  is  based  on  an 
Corrected  Gauss-Sewton  method  combined  with  a  Newton  method  and  an 
Active  Set  method  to  account  for  additional  constraints  on  the 
parameters  to  be  estimated.  The  MLE  approach  offers  the  opportunity  to 
use  a  multi-leg  model,  l.e.  the  platform  to  be  localized  is  assumed  to 
move  according  to  a  piecewise  linear  track,  where  'iach  part  is  referred 
to  as  a  leg.  On  each  leg  uniform  linear  motion  is  assumed.  The  multi-leg 
model  Is  described  as  a  linear  state  model  In  Cartesian  coordinates  with 
non-linear  measurement  equations.  The  MLE  problem,  however,  Is  not 
formulated  in  cartesian  coordinates  since  it  alms  to  estimate  parameters 
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which  are  used  during  real -tine  operation  on  board.  The  method  also 
accounts  for  bearing  and  frequency  neasurement  equations  related  to 
bottom  reflections  of  the  acoustic  signals. 

The  Gauss -Newton  method  requires  only  the  first  order  derivative 
Information  of  the  residuals ,  whereas  the  Newton  method  also  makes  use 
of  the  second  order  derivatives.  In  both  methodes  either  finite 
difference  approximations  or  explicit  analytic  expressions  of  the 
derivatives  can  be  used.  However,  applying  the  explicit  analytic 
expressions  leads  to  a  large  reduction  in  computation  time,  which  is 
much  more  attractive  for  real-time  operation  of  the  method. 

This  report  has  two  objectives.  First  It  recapitalizes  the  theoretical 
background  of  the  MLE  method  In  view  of  the  TMA  problem.  Second, 
analytic  expressions  for  the  first  and  second  order  derivatives  of  the 
residual  with  respect  to  the  target  parameters  are  presented.  The 
expressions  are  required  by  the  Modified  Newton  Method  as  proposed  In 
[Gmellg  Meyling,  1989-1], 

The  report  Is  organized  as  follows .  In  Chapter  2  the  TMA  problem  is 
formulated.  Section  2.2  introduces  the  multi -leg  motion  model,  Section 
2.3  discusses  the  parameters  to  be  estimated  and  their  relation  to  the 
multi -leg  model.  In  Section  2. A  the  measurement  equations  are  formulated 
for  both  direct-path  and  bottom  bounce  propagation. 

In  Chapter  3  the  MLE  method  is  Introduced  and  some  properties  of  the  MLE 
are  discussed  such  as  bias  and  the  Cramer-Rao  lowerbound.  In  Chapter  4 
analytic  expressions  for  the  residual  derivatives  are  determined. 

Readers  who  are  Interested  in  the  results  only  are  recommended  to  skip 
reading  Chapter  4  and  just  to  focus  on  the  tables  of  Chapter  4,  which 
are  recapitulized  in  Appendix  B. 
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2  TMA  MODEL 

2.1  Preliminaries 

Throughout  this  report  the  following  notation  is  used: 

Vectors  and  matrices  are  denoted  In  bold  lower  and  upper  case 
respectively. 

Superscripts  in  parenthesis  denote  leg- Index  numbers. 

Variables  related  to  own  ship  or  target  ship  information  have  the 
superscript  OS  or  TS  respectively, 

-  Similar  superscripts  are  used  for  direct-path  (DP)  and  bottom- 

bounce  (BB)  variables. 

Two  cartesian  coordinate  systems  (X,Y)  will  be  used  where  the  x-  and  y- 
axis  are  related  to  the  geographical  east  and  north.  Absolute  Cartesian 
coordinates  are  defined  such  that  the  origin  represents  a  fixed 
geographical  point,  usually  the  own-ship  position  at  the  time  instant  of 
the  first  measurement.  The  second,  a  relative  Cartesian  coordinate 
system,  moves  along  the  own- ship  track  such  that  the  own- ship  stays  at 
the  origin. 

The  advantage  of  Cartesian  coordinates  is  that  the  absolute  or  relative 
motions  of  both  platforms  can  be  described  straightforward.  So,  the 
relative  uniform  linear  motion  at  time  instant  is  represented  by 
lxk,  yk,  yk]  where  x^,  yk  are  the  relative  position  and  x^,  yk 
are  the  relative  velocity  components .  The  variables  with  superscript  OS 
and  TS  are  absolute  position  and  velocity  components: 

TS  06 

**  -  V  -  V 
yk  -  y?  -  yf 

•  • TS  *06 

\ 

•  •  TS  *06 

yk  -  yk  -  yk 


(2.1) 
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The  position  and  velocity  of  the  target  will  also  be  described  in  polar 
coordinates  [b^,  R^,  Ck,  vk],  where  bk  is  the  absolute  bearing  angle 
taken  from  the  geographical  north  (the  y-axis  of  the  cartesian  system), 
is  the  range  between  target  ship  and  own  ship,  Ck  is  the  target  ship 
course  angle  taken  from  north,  and  vk  is  the  target  ship  speed  at  time 
t^.  Note  that  the  polar  position  coordinates  are  relative  with  respect 
to  the  own  ship  while  the  polar  velocity  components  describe  the 
absolute  velocity  of  the  target  platform.  The  following  relations 
between  Cartesian  and  polar  variables  hold: 


bk  -  AT*H2(xk,yk) 

R*  -  +  y J)1* 

Ck  -  ATAN2(x£s,  yks) 

v*  -  <  +  ykTS’  >H 

where  the  atah2  function  is  defined  by 


ATAN2(x ,  y)  - 


AKTAU  (x/y) 
n  +  ARCTAH(x/y) 
-*r+  A»CTAJi(x/y) 
w/2 
-x/2 


y  >  0 

y  <  0,  x  >  0 
y  <  0,  x  <  0 
y  -  0,  x  >  0 
y  -  0,  x  <  0 


(2.2) 


(2.3) 


Figure  2.1  illustrates  the  definition  of  the  position  and  velocity 
parameters  Just  described. 
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North 


Fig.  2.1:  Position  and  velocity  variables  at  tine  tk 

2.2  Multi-Leg  Model 

The  purpose  of  TMA  is  to  estimate  the  position,  velocity  and  frequency 
parameters  of  an  unknown  platform  from  a  set  measurements  (zk)  at 
measurement  times  (t^.  The  measurements  consist  of  bearing  angles  and 
frequency  measurements  from  the  sonar  system,  which  can  be  obtained  from 
direct  path  or  bottom  bounce  trajectories  to  the  targetship.  Moreover, 
measurements  about  range ,  speed  or  course  of  the  platform  can  be  added 
to  the  set.  In  the  following  we  assume  piecewise  uniform  linear  motion 
of  the  target  ship  (TS)  and  known  position  and  velocity  of  the  own  ship 
(OS). 
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The  TS  notion  can  be  described  by  the  following  difference  equations  in 
cartesian  coordinates: 


where  v^u,  v(1>  are  the  absolute  velocity  components  of  the  target 
ship  which  are  assumed  to  be  constant  during  leg  i.  The  total  number  of 
legs  is  equal  to  m: 


.»  Ju  vk:  5  S,  <  ■  1  -  1 . »  <2  5 

y 

where  r<t‘1>  and  ra)  indicate  the  beginning  and  end  time  of  leg  i,  for 
i  -  1,  . . . ,m.  In  this  report  the  manoeuvre  r(l>  of  the  target  ship  are 
assumed  to  be  known,  i.e.  these  values  may  be  obtained  by  a  manoeuvre 
detection  procedure.  The  time  periods  T^1*  are  defined  as 

Tj1’  -  max(ra)  -  t  ,  0)  +  max(r(t'l)  -  t„,  0) 

K  R  /o  £ 

-  max(r<l)  -  t„,  0)  -  max(r<1_l>  -  tk,  0) 

The  mull- leg  track  is  illustrated  by  Figure  2.2. 
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Fig.  2.2:  The  multi-leg  model. 

The  multi-leg  model  is  determined  by  the  following  discrete-time  state 
equation: 


A  <2-7> 

where 

^  [  xj*  y Is  v'1’  v«\ 


(2.8) 
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« 


N,k 


1 

0 

0 

0 

0 

0 


0 

1 

0 

0 

0 

0 


.(1) 

K, 


0 

1 

0 

0 

0 


0 

T<1) 

»,k 

0 

1 

0 

0 


0 

0 

0 

1 

0 


0 


(2.9) 


0  0 
0  0 


0  0  0  0 
0  0  0  0 


1  0 
0  1 


If  the  model  contains  m  legs  the  state  contains  2(m+l)  variables 
which  have  to  be  estimated.  Additional  to  these  variables  we  also  like 

to  estimate  q  unknown  source  frequency  variables  fQ1,  fw . f0q. 

These  variables  are  assumed  to  result  from  cardinal  peaks  in  the 
frequency  spectrum  of  the  target  platform.  The  corresponding  frequency 
measurements  contain  a  doppler  shift  caused  by  the  radial  motion 
components  of  OS  and  TS.  Section  2.3  discusses  the  corresponding 
measurement  equations  in  detail.  A  frequency  measurement  related  to  time 
tk  and  source  frequency  j  will  be  denoted  by  f  . 

The  state  equation  (2.7)  is  now  expanded  by  q  additional  equations: 


w  -  w  vk  (2-10) 

The  result  is  a  2(m+l)+q  dimensional  state  vector 

*k  ”  I  *k  . v'»>  f01 . f0q  (2.11) 


and  a  transition  matrix 
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I  Ta,.I„  .  ...  Tl“'*I„ 

0 

2  H,k  2  H.k  2 

2  ,q 

0.  X 

0 

2a, 2  2m 

2m,  q 

O 

o 

I 

L  q.2  '  q,2m 

q 

where  I  and  0  indicate  a  n  x  n  identity  matrix  and  a  m  x  n  zero 

n  tn,n  J 

matrix  respectively.  In  the  sequel  of  this  report  we  will  mainly  use  ®k 
which  is  the  inverse  of  k,  i.e.  Tak  -  -  Ta*. 

2.3  Polar  coordinates 

For  operational  usage  it  is  more  convenient  to  use  relative  polar 
position  coordinates  and  absolute  polar  velocity  coordinates  per  leg.  Ue 
introduce  the  polar  state  vector  as 

yk  -  [  bk  ^  C«>  va) .  C‘">  v<»>  f01 . f0?  ]  (2.13) 

with  bk  the  absolute  bearing  angle,  the  range,  Ca>,  v(1>  the  Course 
and  Speed  of  the  platform  at  leg  i  and  fflJ  the  jth  source  frequency. 

The  transition  from  yk  to  yN  can  be  determined  by  using  the  vector 
transformation  function  F:  x  -*  y  and  its  inverse  C:  y  -*  x  : 


yk  -  f(V 
x,  -  G(yk) 


(2.14) 


The  transformations  F  and  C  are  determined  by 
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F  :  bk  -  AiAiu^.y^ 

2  2  U 

\  -  (^  +  y kr 

Cll)  -  ATAH2(v'i,,v‘1)) 

v<i)  _  (v<i>3  +  v'1’2)11 

'  x  y 


i  -  1 . m 


J  -  1 . q 


(2.15) 


G  : 


xk  -  R^sih  bk 
yk  -  Rk  *cos  bk 


r(i) 

X 

-  va)*sin  C(1) 

i  -  1, 

.  ,  m 

r(l) 

y 

-  v(1,*cos  C(1) 

Oj 

1 

Hi 

o 

b*. 

J  -  1.  • 

q 

Note  that  the  own-ship  position  x^,  yk  is  used  implicitely  by  the 
functions  F^)  and  G(yk)  . 

The  transition  from  y^  to  yk  and  vice  versa  is  now  determined  by 


(2  16) 

yk  -  F(*kyG( y„)) 

Note  that  the  own  ship  position  are  assumed  to  be  known  in  order  to  be 
able  to  carry  out  the  transformation. 

2.4  Measurement  Equations 

The  set  of  measured  data  generally  consists  of  bearing  angles  and 
frequency  data  which  result  either  from  a  direct  acoustic  path  (DP)  or 
from  bottom  bounce  reflections  (BB) .  The  acoustic  path  Is  assumed  to  be 
known,  e.g.,  by  using  an  acoustic  propagation  prediction  model. 


TWO  report 


Page 

15 


Moreover,  note  that  the  multi -leg  model  does  not  take  account  of 
variations  in  depth.  Other  type  of  measurements  such  as  range,  course 
and  speed  of  the  target  may  be  exploited  (see  Section  2.4.2).  These 
measurements  may  be  available  from  other  platforms  or  from  information 
about  the  underwater  sound  propagation  conditions.  Section  2.4.1.  deals 
with  typical  passive  sonar  measurements. 


2.4.1  Direct  path  and  bottom-bounce  measurements 

In  case  of  direct  path  acoustics  the  equations  describing  the  measured 

bearings  and  frequencies  are  relative  simple 


2k  “  bk 


+  v 


k 

f  +  v1 

kj  *kj 


(2.17) 


where  bk  is  defined  by  (2.2)  and  f  is  specified  by  the  following 
relation: 


kj 


(2.18) 


The  scalar  c  represents  the  sound  velocity  in  water  (1500m/s)  and  rk 
is  the  range  rate  (the  relative  radial  velocity  component)  at  time  t^. 
The  measurement  noise  is  represented  by  Gaussian  noise  processes  i/£  and 
with  normal  distribution  functions  N(0,ok)  and  N(0,okJ). 

The  second  term  of  the  right-hand  side  in  equation  (2.18)  is  an 
approximation  of  the  doppler  shift.  Here  we  have  assumed  that  the 
measurement  noise  i/kJ  has  a  larger  magnitude  than  the  approximation 
error  of  the  doppler  term.  Moreover,  equation  (2.2)  and  (2.18)  only  hold 
in  cases  where  the  OS  and  TS  do  not  differ  in  depth.  In  case  the  target 
and  the  own  ship  in  reality  move  in  different  horizontal  planes  or  the 
bearing  angles  are  related  to  bottom  bounce  reflections,  the  line  of 
sight  is  projected  onto  a  horizontal  plane.  Figure  2.3  illustrates  the 
underlying  model  for  the  direct-path  and  bottom-bounce  measurements .  The 
elevation  angle  i>  depends  on  the  range  and  the  difference  in  depth  D 
between  TS  and  OS.  The  constant  D  is  either  equal  to  the  real  difference 
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DDr  or  to  the  apparent  difference  in  depth  Dgg  in  case  of  direct  path  or 
bottom  bounce  signals,  respectively: 


D 

D 


DP 

BB 


(2.19) 


where  DDp,  DgJ  and  Dg  are  the  target  ship  depth,  the  own  ship  depth  and 
the  sea  depth.  In  most  long-range  cases  the  difference  DTS  -DQS  is 
negligible.  The  measured  absolute  bearing  angles  bDP  and  bBB  are  obtained 
from: 


l DP  rOS  .  oDP 

\  "  Ck  +  P* 


b“ 


„OS  .  /jBB 
Ck  +  ?k 


(2.20) 


where  /?^p  or  is  the  angle  between  the  projected  line  of  sight  and 
the  OS  course  as  measured  by  the  towed  array  sonar,  fpf  and  /)BB 
depend  on  the  elevation  angle  if.  If  oDP  and  oBB  are  defined  as  the  cosine 
of  the  elevation  angles  ^Dp  and  \fgg  respectively, 


BB  ,  R 

a  -  cos  if  -  — - — c 

BB  (R2+D^  )H 


(2.21) 


then  the  bearing  angles  and  ^B  are  equal  to 


- 


ABCCOS(a^P  COS  Pk)  , 

0  S0k* 

-ABC COS(o£P  COS  0k)  , 

-*  *  < 

«~.ios(<;B  cos  Pk)  , 

0  s  h  * 

-AKCcos(aJB  cos  pk)  , 

-*S  Pk< 

(2.22) 
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where  fj  represents  the  actual  azimuth  of  the  TS  with  respect  to  the  OS 
course : 


bk  "  c?  +  K  (2 . 23) 

In  case  of  large  range  values  the  difference  in  depth  can  be 

neglected.  The  result  is  that  bBP  becomes  equal  to  b^  of  (2.15). 

When  receiving  bottom-bounce  signals  the  frequency  measurements  will 
have  a  different  doppler  component  compared  to  the  direct-path  case.  The 
relative  radial  velocity  between  TS  and  OS  is  represented  by  r.  The 
doppler  shift  is  proportional  to  the  projection  of  r  on  the  bearing 
line  of  the  BB  bearing  angle: 


rDP  -  r  cos  (fi  -  aDP  •  r 


fBB  -  r  cos  ^  -  aBB  •  r 


The  result  is 


kj  Oj 


VI  -  ■«*-> 


(2.24) 

(2.25) 


fkj  -  V1  ~  <V-> 


(2.26) 


where  f  is  the  source  frequency  and  c  is  the  sound  velocity. 
In  case  of  bottom-bounce  the  measurements  are  denoted  by 


<  "  b“  +  < 


kj 


kj 


kj 


(2.27) 


where  i/B  and  w*  are  Gaussian  noise  processes  as  defined  above. 
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Fig.  2.3:  Bottom-bounce  bearing  and  frequency  measurements. 

Note  that  b°p  and  f°p  in  (2.20,  2.22)  and  (2.27)  reduce  to  simple 
range  independent  forms  if  the  difference  in  depth  Ddp  between  OS  and  TS 
is  assumed  to  be  zero.  Moreover,  if  the  angle  between  the  line  of  sight 
and  the  towed  array  is  90  degrees  one  cannot  distinguish  DP  from  BB 
situations . 


2.4.2  Additional  measurements. 

The  additional  measurements  about  range,  course  and  speed  are  denoted  by 


2cu,  _  CU>  +  „c 

zv(1)  -  v(1>  + 


(2.28) 


where  i/J  -  N(0,ej),  -  N(0,o£)  and  -  N(0,ej[)  are  Gaussian 

noise  processes. 
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3  MAXIMUM  LIKELIHOOD  ESTIMATION 

3.1  TMA  formulated  as  an  MLE  problem 

This  Section  deals  with  the  formulation  of  the  general  TMA  estimation 
problem  assuming  Gaussian  distributed  measurement  noise.  Suppose  the 
vector  yH  must  be  estimated  from  the  measurement  set  ZH  which  contains 
all  types  of  measurements  zk  discussed  In  Section  2.4.  Let's  denote  z* 
as  the  variable  to  be  measured  at  tk: 


<  -  w 

zk  “  zk  +  ‘'k  (3.1) 

Z„  -  (zk,  k  -  1 . Ni 

where  i/*  denotes  the  Gaussian  noise  process  N  (0,6*)  corresponding  to 
zk  as  defined  In  Section  2.4.  The  type  of  the  function  hk(yk)  depends  on 
the  measurement  type  at  time  tk,  i.e.  1^  Is  defined  by  (2.20)  through 
(2.23)  in  case  of  a  bearing  measurement  or  (2.26)  in  case  of  a  frequency 
measurement.  The  a  posteriori  probability  density  function  Is  denoted  by 
p(yi|Z)I).  The  maximum  a  posteriori  (MAP)  estimate  Is  found  for  yk  -  yH|R 
such  that  attains  Its  maximum  value.  Using  the  Bayes  rule  and 

taking  the  logarithm  leads  to 

»  P^al2*)  "  «  PC^Iy,)  +  «  P<T«)  ~  «  P(V  (3.2) 

Since  the  last  term  does  not  depend  on  yR  the  MAP  estimate  can  be  found 
by  minimizing  the  function 


i(y„)  -  -l*  p(z,|y,)  -  «  p(y„) 


(3.3) 


TWO  report 


Page 

20 


where  Che  last  tern  represents  the  a  priori  knowledge  on  yg. 

If  there  Is  no  a  priori  knowledge  the  MAP  estimate  reduces  to  the 
maximum  likelihood  estimate,  i.e.  the  last  term  of  (3.3)  is  constant 
when  p(y„)  is  a  uniform  probability  density  function.  Note  that  when  a 
priori  knowledge  is  regarded  as  an  additional  measurement  which  also  has 
Gaussian  statistics  the  MLE  problem  is  in  fact  a  special  case  of  MAP 
estimation. 

The  MLE  problem  is  formulated  as 


min  -ls  p(ZHfyll)  (3.4) 

ys 

Since  the  measurements  in  Zn  are  assumed  to  be  uncorrelated  and 
distributed  according  to  N(0,o*),  i.e. 


P(zjz*)  - 


<2s)'o* 


EXP  (  -H( 


fk^2 


)2  ) 


The  likelihood  function  p(Z|I|y,()  can  be  expressed  as 


(3.5) 


P(ZII|y„)  -  P(zt . z„|y„)  -  £  p(zk|y„) 

...  »  1  *  zv-z*  , 

-  ((2<r)*/Jn  — )  EXP  (-M  2  (-*—*-) a) 

k-l  k-l  £7* 

The  negative  log- likelihood  function  becomes  equal  to 


(3.6) 


V 


+  u  (zwy 


R 

+  £  LR 
k-1 


-w  p(yyp> 


(3.7) 
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where  i(y„)  Is  equal  to 


-  2  (3-8) 

k-1  CTk 

In  the  following  the  weighted  residual  rk(y)  Is  defined  as 


Further  we  denote  the  estimate  of  yk  and  z*  based  on  by  yk|H  and  zkjH 
respectively.  Now  the  MLE  problem  Is  reformulated  as  a  non-linear 
optimization  problem: 


min  (3.10) 

ysls 


A  necessary  condition  for  obtaining  the  ML  estimate  is 


which  leads  to 


(3.11) 


However,  since  the  squared  terms  are  non-linear  In  y,,^  we  have  to  use 
numerical  methods  to  solve  the  MLE  problem. 
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3.2  Numerical  methods 

In  this  section,  only  the  outline  of  the  Gauss-Newton  and  Newton  method 
Is  given.  Extensive  details  can  be  found  in  [Gmelig  Meyling,  1989-1], 

The  Newton  method  uses  the  gradient  S C yH | H >  related  to  some  estimate  yK|M 
to  obtain  a  better  estimate  In  the  sense  that  the  objective  function 
l(Zi,,yI||!()  has  a  lower  value  than  -^^lynia)  • 

The  Newton  method  assumes  a  local  quadratic  behaviour  of  the  objective 
function,  i.e.  the  function  is  written  as  a  Taylor  series  expansion  of 
two  terms: 

i(y  +  s)  =  i(y)  +  g(y)Ts  +  h  sTG(y)  s  (3.12) 

where 


C(y) 


a2i(y) 

3y3yT 


(3.13) 


(See  Appendix  A  for  details  on  vector/matrix  differentials) . 
The  minimum  of  the  right-hand  side  is  obtained  if  s  satisfies 


G(y)*»  -  -g(y) 


(3.14) 


s  is  referred  to  as  the  Newton  direction. 

By  defining  the  N  dimensional  vector  f(y)  containing  the  residuals 
rk(yk|||),  k  -  1.....N  and  the  N  x  n  Jacobian  matrix  J  by 


f(y)  -  trx(y),  r2(y). 


J(y)  - 


at(  y) 
«y 


.  r*<y)]T 


(3.15) 


g(y)  and  G(y)  are  rewritten  as 
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g(y)  -  J(y)Tf(y) 

6(y)  -  J(y)TJ(y)  +  Q<y) 

where  Q(y)  is  equal  Co 


(3.16) 


Q(y)  -  s  r  (y) 

k-1 


32rk(y) 

Syay1 


The  Newton  direction  is  found  by  solving 


(3.17) 


(J(y)TJ(y)+Q(y))*s  -  -J(y)rf(y)  (3.18) 

In  order  to  find  a  descent  direction  the  matrix  G  must  be  positive 
definite.  However,  when  the  residuals  are  large  G  is  not  guaranteed  to 
be  positive  definite  because  of  Q.  A  better  alternative  is  to  neglect  Q. 
Since  JTJ  is  positive  definite,  except  in  special  cases  where  dependency 
between  elements  of  y^  occurs  [Gmelig  Meyling  and  de  Vlieger,  1989],  s 
is  guaranteed  to  be  a  descent  direction.  This  method  is  referred  to  as 
the  Gauss-Newton  method.  A  numerical  robust  way  to  solve  (3.18)  is  to 
use  a  singular  value  decomposition  (SVD)  of  J(y): 

J(y)  -  D  2  V1  (3.19) 

where  0  is  a  N  x  n  orthogonal  matrix,  2  a  n  x  n  diagonal  matrix  and  V  a 
n  x  n  orthonormal  matrix.  The  diagonal  of  2  contains  the  singular  values 
of  matrix  J  In  descending  order: 


2  -  DIAO[  Oj  .  On  ]  , 


2: 


(3.20) 


Moreover  the  following  relations  hold: 
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VTV  -  W1  -  0TU  -  io 

JTJ  -  V  zV  (3.21) 

(JTJ)'X  -  V  Z’2VT  cr4  >  0,  i-1 . n 

The  matrix  V  contains  the  eigenvectors  of  JTJ  and  the  squared  singular 
values  a2  through  o2  are  the  eigenvalues  of  JTJ.  Note  that  If  J  does 
not  have  full  rank  one  or  more  singular  values  will  be  zero  and  a  pseudo 
Inverse  of  JTJ  must  be  used  In  determining  a  search  direction,  l.e.  the 
search  Is  performed  In  a  lower  dimensional  space.  In  fact  the  SVD 
decomposition  accounts  for  the  bad  conditioning  of  the  TMA  problem.  When 
some  singular  values  are  much  lower  than  others  one  can  cancel  those 
directions  corresponding  to  the  small  singular  values  [Gmelig  Meyling, 
1989].  Bad  conditioning  Is  related  to  measurement  noise  and  geometry  of 
the  own  ship  and  the  target  ship  tracks. 

The  outline  of  a  Newton-type  optimization  method  is  as  follows: 

1.  Start  with  an  arbitrary  estimate  y  for  1-1 

2.  Determine  a  search  direction  s  (Newton  or  Gauss-Newton)  by 
solving  equation  (3.18). 

3 .  Perform  a  line  search  along  the  direction  s 

4.  Check  convergence  criteria.  If  not  satisfied  increase  iteration 
counter  1  and  goto  step  2. 

5.  Calculate  the  Inverse  of  G  which  Is  an  estimate  of  the  error 
covariance  of  the  maximum-likelihood  estimate  yj^jj  -  y 

Details  on  Newton- type  methods  can  be  found  in  [Gill,  Hurray  and 
Wright],  Specific  details  on  Newton-type  methods  for  TMA  problems  are 
described  in  [Gmelig  Meyling] . 
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3.3  Unbiased  Estimates  and  Che  Cramer-Rao  lower  bound 

The  estimation  error  e(Z)  is  defined  by 


•<z>  “  r, I,  -  y„  0.22) 

Any  unbiased  estimation  y^  of  the  parameter  vector  yN  is  characterized 
by 

E(e(Z)|yH)  -  Jze(Z)  p(Z|y„)dZ  -  0  (J23) 

where  E(e|y)  denotes  the  mathematical  expectation  of  e  knowing  y.  If  we 
consider  yR|S  as  an  biased  estimate  of  yH  with  bias  A(yN)  the  expected 
value  can  be  written  as 


E(e(Z)|y„)  -  A(yn) 

Differentiation  with  respect  to  yM  leads  to 


3  u(p(Z|y,,)) 


P(Z|y„)dZ 


I  + 


n 


3X(y„) 


Now  the  error  covariance  matrix  At  is  defined  as  follows 


(3.24) 


(3.25) 


\  - 

-  El[e(Z) 


x(y,)He(Z)  -  A(y„)lIp(Ziys)dz 
My„)ll«(Z)  -  A<y,,) ]T|y„> 


(3.26) 


The  vector  3ui(p(Z|y||) )/5y1|  is  a  stochastic  vector  with  expectation  equal 
to  zero  and  a  covariance  matrix  equal  to  the  matrix  K: 
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(3.27) 


The  n  x  n  matrix  M  is  called  the  Fisher  Information  matrix.  The  Cramer- 
Rao  theorem  [v. Trees,  Eykhoffj  states  that  for  any  unbiased  estimate  of 
yNjN  the  following  inequality  holds  (see  appendix  C)  : 

A#  2:  M'1  (3.28) 

The  meaning  of  this  theorem  is  that  for  each  unbiased  estimate  the  error 
covariance  matrix  has  a  lower  bound,  which  is  the  inverse  of  the  Fisher 
information  matrix.  The  lower  bound  can  be  determined  exactly  if  the 
real  target  state  ys  is  known.  When  Monte  Carlo  experiments  are  used  to 
show  the  performance  of  the  MLE  method  it  is  thus  possible  to  determine 
the  lower  bound.  When  the  estimation  method  is  efficient  the  error 
covariance  matrix  is  equal  to  the  lower  bound.  The  error  covariance 
matrix  can  be  either  estimated  from  the  Monte  Carlo  simulation  results 
or  by  using  the  matrix  (J(yH|H)Tj(y1I|„)+<)(y(I|K))'1-  The  estimated  error 
covariance  matrix  can  be  compared  to  the  lower  bound  in  order  to  get  an 
efficiency  measure  of  the  MLE. 

Using  the  definitions  of  g  and  G  the  matrix  M  can  be  rewritten  as 

H  -  E{g<yK)g(yIt)1  |yH>  -  Jiy./JCy*)  (3.29) 

Unfortunately  it  cannot  be  proven  that  y^  is  unbiased.  Moreover  one  is 
not  able  to  determine  an  analitlc  expression  of  the  bias  My„)  ■ 
Experimental  results  in  [Gmelig  Meyllng  and  de  Vlieger,  1989-2]  show 
however  that  in  many  cases  the  MLE  will  be  unbiased  and  one  can  use  the 
Cramer-Rao  lowerbound  to  verify  the  efficiency  of  the  MLE  method. 
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3.4  Quality  of  Fit  and  Confidence  regions. 

Given  an  estimate  y  the  likelihood  was  defined  by  (3.6)  and  (3.8): 


p(Z|y)  -  (  (2Jt)'n/2  H  (oVl>  BP(-hi(Z,y)) 

k-1  ' 

.  (3.30) 

N  Z  -Z„ 

i(Z  y  )  -  S  (-*—*)’ 

"  k-i  o' 

The  quantity  i(Z,y)  at  its  minimum  is  distributed  acc<  ding  to  a  chi- 
square  distribution  function  for  N-n  degrees  of  freedom,  where  n  is 
equal  to  the  dimension  of  y.  Suppose  we  have  a  realization  ^  (a  set  of 
measurements)  and  an  estimate  yMLE.  Other  possible  realizations  Z  which 
lead  to  other  MLE  estimates  y  with  a  likehood  i(Z,y)  such  that 

i(Z,y)  >  i(ZR,yMLI)  (3.31) 

result  in  MLE  solutions  which  do  not  fit  the  data  Z  as  well  as  yMLE  fits 
ZR.  Hence  the  probability  Pr  ( i(Z,y)>i(ZR,yMLE)  can  be  used  as  a 
quantitive  measure  for  the  quality  of  fit  of  yf1Lt: 

Pr(i(Z,y)>i(ZR,yMLl))  -  Q(J^B  ,  h* (Z^y^) )  (3.32) 

Q(v,a)  is  refered  to  as  the  incomplete  gamma  function: 


Q(v,a)  -  -1-  f  xv  l  nc?(-x)dx  (3.33) 

r(v)  Ja 

If  the  measurement  noise  realizations  resulted  into  realization  Z^  with 
a  low  probability  of  occurrence  the  likelihood  of  the  ML  estimate  will 
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be  low  and  the  probability  1-Pr  will  Indicate  that  there  Is  a  high 
chance  chat  someother  realization  would  fit  better.  However,  the  quality 
of  fit  may  not  be  confused  with  the  confidence  in  the  estimate.  One  may 
obtain  a  bad  fit,  for  instance  because  of  a  bad  measurement,  but  still 
obtain  an  excellent  estimate  and  visa  versa. 


Let's  now  consider  all  possible  estimates  y  for  which  the  likelihood  is 
larger  or  equal  to  a  lowerbound  a  p(Z-ZR|yMLI) ,  for  0  5  a  s  1: 

a  p(Z-ZR|yMLE)  <  pCZ-ZJy)  <  pCZ-ZJy"1-1)  (3.34) 

Hence  we  select  all  estimates  around  for  which 

f<zR.yMLI)  s  fCZR.y)  s  f(ZR,yMLE)  -  2u«  (3.35) 

The  set  (y)  satisfying  condition  (3.35)  defines  a  likelihood  region  R  in 
the  state  space.  The  probability  that  y  lies  in  the  region  R  given  the 
measurement  set  Z^  is: 

PrlyeSIZ-ZR)  p(y|Z-ZR)dy  (3.36) 

By  using  the  Bayes  rule  we  obtain 

p^z-z*)  -  _g<z-AIZ>?^ 

^  ;  p(Z-ZR|y)p(y)dy  (3.37) 

y 


Note  that  the  nominator  of  (3.37)  is  a  normalization  constant.  The 
desired  probability  can  be  calculated  if  p(y)  is  available.  Remember 
that  p(y)  represents  the  a  priory  knowledge  about  the  target  track  and 
that  the  KLE  method  does  not  account  for  this  knowledge .  The  worst  case 
situation  is  at  best  represented  by  a  uniform  distribution: 


p(y) 


1/A”  -A/2  5  yt  s  A/2,  i-l,..,n 

0  yt  S  -A/2  or  S  A/2 


(3.38) 
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Hence,  for  A-*®  the  a  priori  knowledge  vanishes.  Substituting  (3.38)  Into 
(3.37)  results  In  a  probability 


PrCyeSIZ-Zj^) 


Jgp(z~zRly)dy 

jyP(Z-ZR|y)dy 


(3.39) 


Substitution  of  (3.30)  leads  to 

-  hi  (ZR,y)  )dy 

PrlyeSIZ-2^)  -  (3.40) 

J  otP( -hi(ZR,y)  )dy 


which  is  the  probability  that  the  real  target  parameters  lie  in  a 
confidence  region  R  defined  by  (3.35)  based  on  the  knowledge  of  the 
measurement  set  ZR.  In  general  the  boundary  of  R  will  not  be  an 
ellipsoid  because  of  the  nonlinear  behaviour  of  the  likelihood  function. 
The  shape  of  R  can  be  found  by  using  a  Monte  Carlo  Integration  approach, 
i.e.  by  producing  a  sufficiently  large  number  of  random  vectors  y,  which 
are  drawn  from  a  uniform  distribution  function  that  fully  covers  the 
region  R,  one  can  display  those  vectors  y  graphically  that  satisfie  the 
condition  (3.35).  Simultaneously  one  can  approximate  the  Integrals  of 
(3.40)  numerically  by  using  both  the  accepted  and  the  rejected  vectors. 
The  displayed  points  are  scattered  in  the  region  R  and  hence  show  the 
shape  of  R. 


4  RESIDUAL  DERIVATIVES 

The  application  of  Newton- type  algorithms  implies  the  usage  of  first  and 
second  order  derivatives  of  the  estimate  with  respect  to  the  state 
estimate  y„|„.  This  Chapter  is  concerned  with  the  determination  of  these 
derivatives.  First  we  focus  on  bearing  and  frequency  residuals  with 
respect  to  direct-path  and  bottom-bounce  measurements.  Then  range, 
course  and  speed  residuals  are  considered  in  Section  4.2.  In  Section  4.3 
second  derivatives  of  the  residuals  are  determined. 

4.1  First  derivatives  of  bearing  and  frequency  residuals 

The  general  vector  form  of  the  first  derivative  of  a  residual  tk(yk|„)  is 
equal  to 


3rk^ykllP  _  _  1  *Vr<*k.^<W ) )  ^  -y 

3yN|H  ayn|H 

Throughout  this  Chapter  we  just  concentrate  on  the  derivatives  of  hj, 
instead  of  rk(yk|H)  •  According  to  the  rules  of  Appendix  A  the  following 
expression  is  obtained: 


ahklZn->  -  ac(W-t»  .  af<W,  ahk(ykiw) 

3yS|»  3y*|»  3*k|*  3yk|* 

which  can  also  be  denoted  as 


dZll9  _  .gi  .Km 

3y»|»  8yS|»  ’  3xk\» 


where  z*,,  -  hk(yk|>) . 
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In  most  cases  computational  efficiency  Is  Improved  by  calculating  the 
derivatives  from  explicit  analytic  expressions  Instead  of  using  (4.2)  or 
(4.3).  However  (4.2)  and  (4.3)  will  serve  as  a  guide  to  derive  the 
individual  expressions  of  each  element  of  each  type  of  residual 
derivative.  First,  consider  the  direct-path  bearing  estimate  bk|K.  The 
result  of  equation  (4.3)  is 


abtl»  _  a>KlH  ,t  3«AS2(xi|M,yk|M) 
k'K  a*k|H 

Further  we  have 


Ruin035  bs|s  -KhIs81"  bs|» 
s»  bH(H  COS  b#,„ 


3*S|S 

3y«|H 


°2,Z 


'  V(1)C0S  C»)  -v(i)si*  CU> 
SIS  c<t>  cos  cu> 


0 


<J.3 


\ 


(4.5) 


Note  that  (4.5)  Is  a  block  diagonal  matrix  with  m+1  subaatrlces  of 
dimension  2x2  and  one  q  x  q  Identity  matrix  on  its  diagonal.  Therefore 
it  is  much  more  convenient  to  rewrite  (4.4)  as  a  number  of  expressions 
with  respect  to  bearing,  range,  speed  and  course  of  each  leg,  and  each 
source  frequency: 
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The  eetioated  direct-path  frequency  equation  is  denoted  by 
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kj|* 


(A. 7) 


where  ^0j|«  is  tbe  esC^mate  of  source  frequency  f  as  a  result  of  the 
set  2^.  The  derivative  with  respect  to  y^  can  be  derived  according  to 
(4.3).  It  is  however  easier  to  obtain  the  derivative  straight-forward 
from 


£fj. 


ilia _ klr  +  8tm[i»  (i  _ 


af„ 


S|H 


c  a  j. 


K|B 


ay 


(4.8) 


si* 


Note  that  is  equal  to  a2(otl)+J  where  et  is  the  identity  vector 

with  a  unity  element  at  the  itb  position.  Substitution  in  (4.8)  leads  to 


8ffcll»  _  _  _!si  +  (1  _  iklS.)  , 


2Cb+1)+J 


ay.i*  c  ay«is 

The  first  term  of  (4.9)  depends  on  the  range-rate  derivative  rk|N: 


(4.9) 


*k|«  -  <v“’  -  ^>SII,bs|»  +  <<’  -  >«»  bk|, 


tu'l)  s  tjj  <  ttl)  or  tk  -  t(“’ 


For  the  sake  of  convenient  notation  we  introduce  vjj^J  as 


vi?»  ”  (vi1J  "  *?>«*  bk|,  ~  <v"'  ~  y">«»  \ 


(i)  _  AOSv 


t(1-»>  s  t*  <  t(1>  or  -  t‘"> 


(4.10) 


(4.11) 


Equivalent  expressions  are  : 
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*1* 


k|» 


Ci|i  -  *?)«*  bHH 

+  ciii  -  yf)cos  bM« 

KJ>  -  *?)«•  bkl» 

-  -  y>"  bk!» 


Differentiation  of  (4.10)  with  respect  to  y^,,  leads  to 


3*kls 

TA*  at>k|H 

vm 

fills 

h)h' 

ac‘f> 

g  klK 

ay*i* 

k|"  ay«in 

Vk|S 

SlH(bk|H  - 

ay«iK 

+  cos(bk|ll  - 

fills 

3y*is 

av™ 

^hin 

kl" 

v(l> 

k|s 

““(hklm  ~ 

fills 

3y»is 

-  SI*(bklH  ~ 

fills 

c*l»> 

av(1) 

!lsls 

3y*n. 

Note  that  C‘|’)  -  C*j’)  £«  «H  k  since  the  TS  course  at  leg  i 
remains  constant.  The  result  is  a  function  of  the  bearing  derivatives: 


-  *S  -  <4p 

dy*|»  "y*!* 

+  costb  ,  —  C<‘>)  e 

“"'“kl*  »|s'  21+2 


3v 


rTAK 

VlP 

kl"  3y.i. 

+  vkl* 

-  s«<bk|s  - 

C*|*'  *2 

b  -  f<‘h  a 
°k|*  C*|»'  *21+1 


(4.13) 


Substituting  the  bearing  derivatives  of  Table  4.1  in  (4.13)  leads  to  the 
result  summarized  in  Table  4.2. 
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4.2  First  derivatives  of  range,  course  and  speed. 

The  first  derivatives  of  the  course  and  speed  estimates  are  obtained 
straight  forward: 


ac 


ay, 

dv, 


If 

(1) 


(4.14) 


*7. 


•  I* 


•  th 


where  Cj  Is  the  Identity  vector  with  a  unity  element  at  the  j 
position. 

The  Range  derivatives  are  obtained  by  applying  the  general  equation 
(4.3): 


k-N  a*kis 


(4.15) 


where 


-  t  sin  b  I  cos  b  .  0  .  0  J1 


(4.16) 


Carrying  out  the  matrix  multiplications  leads  to  the  results  in  Table 
4.3. 


The  variable  is  defined  as 
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in  correspondence  to  the  definitions  in  (2.21)  and  (2.19)11. 
Similar,  the  bottom-bounce  frequency  estimate  is  equal  to 


kj 


-  f 


Oils 


(1  - 


lilih 


(4.19) 


The  derivative  of  bBB  can  be  expressed  as 


3bBBK  3akccos(u)j- 

3yB|N  8u 


r  w1"  0*i» 


3K 


ay, 


»|H 


8aklH  3RklS 


3  R, 


k|M 


] 

(4.20) 


where  u  -  cos  /9kjK  and 

3  /jiccos(u)  -1 

3u  (l-oJ|8.cos2(bk|K))l'J 

Further 


(4.21) 


da. 


_B£_ 


3Rk|H  (RCis  + 


dsb> 


3/2 


(4.22) 


The  range  derivative  in  (4.20)  can  be  found  from  table  4.3.  Hence  the 
result  can  be  found  from 


3y*i» 


ab“„  abn„  +  aAis 

abk|»  a3rk|m  a“k|»  aRk|»  **k|« 


(4.23) 


1)  Rwiwfctr  that  m  still  uium  D  to  ba  taro.  Not*  that  whan  la  substantial 

TS  08  TS  06 

larga,  wa  alao  hava  to  rafonwlata  tha  dlract-path  bsarlng  and  traquancy  aquations  by 
ualnc  tha  botto*'bounea  aquations  and  raplaca  Dgg  by  D^p. 


where 


3b; 


n 


3b  , 
“kin 

Him 

Him 

Him 


Wx»Him>-51<”Hh> 

(KiXfbkl,))1'2 


^HlM^2 


(iHiM-^Hi.)) 


1/2 


H2|M  +  D~> 


fo312 


(4.24) 


The  direct-path  bearing  and  the  range  derivatives  can  be  found  from  the 
tables  4.1  and  4.3  respectively.  There  is  one  trivial  case: 


-  0 


3f 


j  -  1.  ....  q 


(4.25) 


OJ  |H 
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fkj|H  ”  f0j|H  +  “klN^kjlK  "  ( A  .  26 ) 

Remember  that  fkJ|N  is  equal  to  if  DDp  is  zero! 

The  first  derivative  can  be  calculated  straightforward  by  using  Table 
4 . 2  and  4.3: 

3fk.)lH  _  0  3fk.tlH  _  f0jlHrklS  3°k IW  3RklW 

ayK|N  k|N  ^NlH  c  3Rk|N  ^KlH 

af  (427> 

+  (l^klB)^i» 
klH 


The  result  is  shown  in  Table  4.5. 


Table  4.5  First  Derivatives  of  Bottom- bounce  frequency. 
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4.4  Second  derivatives  of  range,  bearing  and  frequency  residuals 

The  second  derivatives  of  direct  path  and  bottom-bounce  bearing  and 
frequency  estimates  are  obtained  from  the  Tables  4.1  through  4.S  of  the 
previous  sections.  First  we  focus  on  direct  path  derivatives.  From  the 
first  order  derivative  5bk |R/3bR||(  in  Table  4.1  we  obtain 


3zb 


iiN_ 


_  SslB 


.3RHh 


d-db 


N  |  N 


R2  C0S<bH|H-  \\t>)-rLS  +  T~  C°S<bll|K- 
*|h  a 


,glJkl»  _  ^gJJ! 


d  • 


+  SXH(b,||H-  b..H)(— ^  ) 

"|li  k|H  a •  a • 


(4.29) 


where  ♦  indicates  one  of  the  estimated  variables  in  y„|N-  Note  that  the 
second  term  is  only  non-zero  if  •  indicates  R„|H-  Similar,  when  •  is 
equal  to  bH((i  the  term  with  db^/3*  will  become  non-zero.  Many  of  the 
terms  can  be  expressed  easily  by  using  first  order  derivatives.  In  this 
way  the  second  derivatives  can  be  expressed  by: 


a\lH  .  _  fill  "klf  +  1  3bkls 

a*  abM|H  ^IK  abH|K  a*  ^1*  a* 


+ 


^IS 


ab>i»  (3bn„  _  abwm  } 

aR»l»  a’  a’ 


(4.30) 


This  procedure  is  repeated  to  obtain  second  derivatives  from  all  the 
first  derivatives  of  Table  4.1: 
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a2b 


klH 


1 

Kt* 


aRVlN 

2  bk|.)^7 


1  3b  , 

-^.C0S(b»,«- v)(i^* 


£5kn 


3zb, 


ilB_ 


w(i)  an 

-  -  5s  VN|i  “(W  S|i> 


K\* 


rU> 


^Is 

'r(i) 

h 

S|. 


b|m 


sm(b 


0< 


k|N  *|N 


3b, 


3vif: 

uk|M~  a. 


+  -U  cos(b  ,  -  C‘t>)  ->«> 


8C  I 

_  °_NJ.K  ) 

8* 


3zb, 


iia_ 


T*(i  ) 

-fcJL  «»<W  C,‘l> 


Si. 


k|H  *1.'  3. 


)3AiH 


T‘“  ...  3b. ,  3C“> 

- ***  cos(b  |  -  C‘{’)(— ^ - ***  ) 

k|H  *l"  3*  3* 


3-3f, 


ill - o 


oj 


(4.31) 


By  using  the  first  order  derivatives  expressions  as  much  as  possible  the 
result  in  Table  4.6  is  obtained. 
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Table  4.6  Second  Derivatives  of  Direct-Path  bearing. 

In  the  following  the  second  order  derivatives  of  the  direct-path 
frequency  estimate  fkJ  are  obtained  In  a  similar  way: 
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a2  f. 


k.llN 


a  f. 


af, 


ktiN  g '•!!.<  |H  +  3fkllN  PVH[K  _  '•Qllll 


3V»> 


3- 


f0j|.  3CK|i  3 


3btlu  3v; 


K\l  K\l  3* 


C  3CJ|J  3* 


MiS. 


,  yli)  3fHlll  /3bkll»  _  3CSlll  ■.  +  f 0.1  la  Vklu  3bklH  3RklW 

-  »*  3.  3. 


7(1) 

'*1* 


3v 


(1) 


Si.  3Ciii3' 


(A. 35) 


Similar,  a^f^iH/ayuijjavJj’  is  stated  below: 


3zf, 


idle _ _ 

o__( 1 ) 

3v*|a 


af, 


ma 


ojls 


3v»|a 


3  f  o.i  In 


a  c, 


HlN 


+  f07lH  VklH  3b>lN  3Rkl„ 

c  S|a  3vJ|h  3* 


foi  Iw 
c 


av™1 

kls 


(4.36) 


The  derivatives  3zfkJ |H/ayH|„3 f0J|i,  are  obtained  straight  forward  from  the 
range-rate  derivatives  (4.13).  The  result  is  summarized  in  Table  4.7. 
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Table  4.8  Second  order  Renge  derivatives. 

The  bottom-bounce  bearing  derivatives  are  obtained  by  differentiation  of 
the  expression  (4.23): 


3* b: 


klB_ 


8-  8y. 


»|» 


+  J«a_  +  3J^Ll . ^J1L  aAin_ 

3*  abki»  3*«i«  abki»  3*33rsi»  3'3aki*  3Rki»  3y»i» 


,  3b”,  3Sl,  3RH„  ,  3b”„  3°H,  3\l, 

3aki*  3#3Rki»  37*i»  3ai,i»  3Rkis  3*  3y»is 


We  memorize  that  the  first  order  derivatives  in  (4.38)  are  already  known 
from  the  Tables  4.1,  4.3  and  4.4.  The  derivatives  3ab®J|/3»flb^|1|  and 
a2b"./3*3«*k|.  can  be  found  from: 
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TTaT  "  (7T  (1'^i»“s  >aki«sl,^ki.)SIO'<^i«) 

k|» 


+  <Hj|I/W,^|B>'1/#  —  (“klH  ««(^lB)) 


(4.39) 


a*  a^,,  a- 


<17*  (1-‘Sc|»C0S  ^k|n>*  )-««C/3k|ir> 


Ic  Is  obvious  chat 


—  <1-^2  cos^  )-"2 - y  ^ 

a-  ^l*  k!"  a^^cos2^,,/3'2  a- 


sii,4ih_  8Aih 


(4.40) 


(4.41) 
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SubsCltucion  in  (4.39)  leads  Co 


aZbki*  sig«(A  ,  )  2  2 

*<P&* 1  {^-^- 


+  <1^i*co8Xi*)si''(^ki*))jri!L  +  <-^t*cos^kiH>s«2<^ki»> 


+  (l-ok„,cos2ak|l,)<*k|Hco8(^|1|))-^  ( 


(4.42) 


SI0*(^kl»)  0akl*  5  0t5kll 

“  (i-m,^,,)*/*  ,SI’(/Wa.  +  “ki^kisd-M*)^ 
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rS  -  <1^^Xi»>'3,218i"<*m.>j^  -  “k|«,cos3^k|»)rrL")  <*•«> 


3«  3oj 


The  2Bl1  derivative  |„/a •  3R^ (l|  in  the  4th  term  of  (4.38)  is  equal  to 


JSii.  -3D»  Si»  ats 

fl-3Rk|s  <«J|.  +  d»b)S/2  a‘ 


The  result  is  recapitulized  in  Table  4.9. 


(4.44) 


a2bkf*  _  >  3_$ls 

a-  abkis  <Ki^i.)3/2  kl"  a«k|s  a* 

3b 

+  °k[sC0S^(«(1-°k2|s)-^11!) 

57^,;  (1-^ihcos2^i.>'3/2<si"^i.>^- 

aSls  -3Psb  »n,  aRkls 

a‘aRk|s  <*J|.  +  I4)S/I«- 

a2b“s  _  a2b“,  abH,  ,  abklB  a\lS  ,  aX?H  a°klH  aRklS 

a-  ay^  a-  abk|„  ay,,,  abk|g  a>aym  a*a«^|g  a**,,,  ayg|K 

,  gSl  3Rk|>  ,  ab“k  3W*  a\lK 

a<si»  a*  3Rkis  ay«is  aaki»  aRki,  a-  ay»i» 


L 


Table  4.9  Second  order  bottom-bounce  bearing  derivatives. 

The  bottom-bounce  frequency  derivatives  are  determined  from  (4.26). 
By  partial  differentiation  of  fJJ|g  with  respect  to  f^,  fkj(l(  and 
we  obtain 
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(  iSiliL  _  d f o.i In  ^“klM  _  f  *xln  a2°fc|K 
fl.  3.  a^,  0JlH  c  a.  a^. 


(4.47) 


Note  that  terms  with  afW|S  /  ^  become  zero  since  we  use  the  partial 

derivatives  of  (4.26)  with  respect  to  fkjl„  and  to  obtain  (4.45). 
The  result  is  stated  in  Table  4.10. 
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5  CONCLUSIONS 

In  this  report  the  outline  of  the  MLE  method  for  TMA  is  described.  The 
method  requires  some  kind  of  optimization  procedure.  Since  TMA  problems 
are  often  ill-conditioned  a  robust  numeric  optimization  method  is 
required.  In  [Gmelig  Meyling]  a  Newton-type  method  is  proposed  by  which 
TMA  problems  can  be  solved  very  efficiently.  This  report  describes  how 
numerical  experiments  can  be  checked  by  using  the  Cramer-Rao  bound  and 
moreover  that  the  MLE  may  not  produce  unbiased  estimates.  Newton-type 
methods  require  first  order  derivatives  of  the  function  to  be  optimized. 
Near  the  optimum  also  second  order  derivative  information  can  be  used  to 
improve  the  convergence  of  the  method.  Although  derivative  information 
may  be  obtained  by  using  finite  differences  much  computation  time  is 
saved  by  using  the  analytic  expressions  which  have  been  derived  in 
chapter  4.  The  experimental  results  and  specific  details  about  the 
Newton- type  method  which  we  prefer  can  be  found  in  two  related  reports 
[Gmelig  Meyling,  Gmelig  Meyling  and  de  Vlieger] . 
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SOME  MATRIX  DIFFERENTIATION  RULES 


Suppose  x,  y,  s  and  m  are  arbitrary  column  vectors,  M  is  a  matrix,  and  f 
is  a  scalar  function,  then  the  following  theorems  hold: 

(1)  If  y  -  Mx  then 

dy  -  M dx  implies  dy/dx  -  M'  and  conversely 

(2)  If  f(x)  -  m'x  then 

d f  -  dm'x  -  m’dx  implies  df/dx  -  m  and  conversely 

(3)  If  x  -  mf  then 

dx  -  mdf  Implies  dx/df  -  m'  and  conversely 

(4)  dz/dx  -  dy/dx  .  dz/dy 

(5)  If  f(x)  -  x’Mx  then 

df  -  (dx) 'Mx  +  x’Mdx  implies  df/dx  -  (M+M’ )x 
and  conversely. 

(6)  df/dxdy  -  d(df/dy)/dx 

(7)  If  f  -  x'My  then 

df  -  (dx)'My  +  x'Mdy  Implies 
df/dy  -  M’x  and  d*f/dxdy  -  M 

(8)  If  f  -  x'Mx  then 

d’f/dx2  -  M'+  M 
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TABLES  OF  CHAPTER  4 


In  this  appendix  the  tables  with  first  and  second  order  derivatives  are 
repeated  as  a  quick  look  facility  for  those  readers  who  are  only 
interested  in  the  results  of  Chapter  4. 


Table  4.1  First  derivatives  of  Direct-path  bearing  estimate  b..„. 
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Table  4.2 


First  derivatives  of  dlrect~path  frequency  estimate  fkJ|B. 
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Table  4.3  First  derivatives  of  Range  estimate  R ^ 
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Table  4.4  First  Derivatives  of  Bottom-bounce  bearing. 
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Table  4.5  First  Derivatives  of  Bottom-bounce  frequency. 
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Table  4.6  Second  Derivatives  of  Direct-Path  bearing. 
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Table  4.7  Second  Derivatives  of  Direct-Path  frequency 
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Table  A. 8  Second  order  Range  derivatives. 
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Table  4.9  Second  order  bottom-bounce  bearing  derivatives. 
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Table  4.10  Second  order  bottom-bounce  frequency  derivatives. 
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This  theorea  can  be  proven  by  defining  the  vector  X  as 


X 


•  (Z)  -  A(y„) 


dm(p(Z|y„)) 


The  covariance  matrix  of  X  must  be  non-negative  definite: 


(C.l) 


e ( xx1 1 y„ ) 


I  + 


n 


dACy,)' 


M 


a  0 


(C.2) 


Hence,  for  any  vectors  x  and  y  the  following  quadratic  fora  must  be 
nonnegative : 


T  t  3My„)  ,  9X(y.)  . 

xTA  X  +  XT(I  +  - L-)y  +  yt(I  +  - L)IX  +  yTMy  2:  0  (C.3) 

3y„  dy. 

Assuming  that  both  A(  and  M  are  positive  definite  the  minimum  of  the 
left  expression  Is  attained  for: 


-l  3A(y  ) 

y  -  -  *  "a,  ♦  -r-*-)  x 


(C.4) 


Substituting  (C.4)  into  (C.3)  leads  to  the  following  condition: 
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xT[  A.  -  <!„ 


ax(  y,) 


ax(y.) 


ay. 


)T  ]*  *  o 


(C.5) 


Hence,  when  the  estimate  has  a  bias  My„)  the  Cramer-Rao  theorem  changes 
Into: 


A.  *  + 


3A(y[)) 

ay* 


n 


+ 


iiiZal)  t 
ay» 


If  A(yK)  is  zero,  condition  (C.6) 


reduces  to  (3.28).» 


(C.6) 
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